This effort is also part of the VAIA FRONT project - FRom lessong learned to future Options .
We have to assign land-cover attributes to nodes over a large area for simulating the impact of strong winds over forest areas.
Several layers are processed to aggregate information over nodes of a regular grid.
Data were processed with Google Earth Engine (GEE) and R.
Code is reported below (snippets) and available in Github .
Outputs are available at Github page in the folder “out”.
Training and validation was done over areas that were defined over areas that were damaged from the VAIA storm event.
The data is extracted from the windthrow database that was created for the paper by (Forzieri et al. 2020) linking to an open database containing polygons of damaged areas. It must be noted that some damaged areas were left out because only slightly damaged, or with small areas or outright missed during the visual interpretation process. Therefore we should be careful to conclude that all areas not covered by the polygons are un-damaged. Nevertheless these data provide us with samples of damaged areas.
There are 7416 polygons in the dataset: 5235 have information on the degree of damage, expressed between 0 and 1, 2181 have NULL value (expressed as -9999).
A regular grid with 1 km spacing was used over the study area, for a total of 97524 nodes.
N.B. A nice alternative could be to have an hexagonal grid, which has several geometric advantages.
Cannot convert our regular grid to hexagonal grid because hexagonal requires offset at alternate rows.
Grid is obviously not aligned to a regular latitude longitude grid as distance between nodes would decrease moving away from the equator.
A custom projection Lambert Conformal Conical projection was designed (Marika Koukoula) with the following PROJ
myproj<-"+proj=lcc +lat_1=45.827 +lat_2=45.827 +lat_0=45.827 +lon_0=11.625 +x_0=4000000 +y_0=2800000 +ellps=GRS80"
Each node was expanded to a square
st_bbox_by_feature = function(x) {
x = st_geometry(x)
f <- function(y) st_as_sfc(st_bbox(y,), crs=myproj)
do.call("c", lapply(x, f))
}
nodes.myproj.buffered<-st_buffer(nodes.myproj, 500, nQuadSegs = 1)
squares<-st_bbox_by_feature(nodes.myproj2)
squares <- squares %>% st_set_crs( myproj)
The square lattice was used for further processing with Google Earth Engine and with R.
The 97524 square polygons were uploaded to Google Earth Engine for processing raster products derived from satellite data.
The grid was used to extract of the following two files:
copernicus.csv
Dictionary with histogram of class frequency of Copernicus classes inside the 1 km x 1 km area square around the nodes of the grid (domain?). Extracted using Google Earth Engine map/reduce over 2018 COPERNICUS land cover
nodesWithGEEvars_crsXXXX
XXXX = 4326 (latitude longitude in WGS88) or XXXX=LCCcustom (custom Lambert Conformal Conical see here ) ) with percentiles [10,25,50,75,90] of the following as attributes. A shapefile with nodes and the five percentiles for each of the following attributes.
Tree canopy cover for year 2000, defined as canopy closure for all vegetation taller than 5m in height. (https://developers.google.com/earth-engine/datasets/catalog/UMD_hansen_global_forest_change_2018_v1_6) Hansen, Potapov, Moore, Hancher et al. “High-resolution global maps of 21st-century forest cover change.” Science 342.6160 (2013): 850-853.
Tree canopy height for year 2005, defined as canopy closure for all vegetation taller than 5m in height. (https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2011JG001708)
Aspect
Slope
DEM (Height)
The last three from the SRTM mission. Slope and aspect in degrees. Height is meters A.S.L. (https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2011JG001708)
The following code was used in GEE.
var getCentroid = function(feature){
return feature.centroid();
};
var slp = ee.Terrain.slope(srtm).clip(domain_buffered.geometry().bounds()) ;
var asp = ee.Terrain.aspect(srtm).clip(domain_buffered.geometry().bounds()) ;
var dem = srtm.clip(domain_buffered.geometry().bounds()) ;
var image = hansen.select(['treecover2000']).rename(['cCov'])
.addBands(canopyHeight.rename(['cHgt']))
.addBands(slp.rename(['slp']))
.addBands(asp.rename(['asp']))
.addBands(dem.rename(['dem'])) ;
var treeCover = image.reduceRegions({
reducer: ee.Reducer.percentile([10,25,50,75,90]),
collection:domain_buffered
});
var unBuffer = ee.FeatureCollection(treeCover.map(getCentroid));
Export.table.toDrive({
collection: unBuffer,
description:'mapRedVars',
fileFormat: 'SHP'
});
Processing 97524 squares intersecting them with 121924 polygons of forest category and 7416 polygons of damage areas too 120 seconds with 14 CPU threads parallel tasking (not sure if this is the correct terminology - not a parallel guru).
Table with variables: asp=aspect, slp=slope, dem=DEM, cHgt=canopy height, cCv=canopy cover, mainCat=main category, mnCt_Cv=main category cover, dmg_Cov=damage cover, dmg_wCov=damage weighted cover
All cover values are in m^2 and can be converted to relative cover dividing by the square area (1 km^2).
| Column Name | Type | Min | Mean | Max |
|---|---|---|---|---|
| asp_p10 | numeric | 0 | 71.05 | 336.79 |
| asp_p25 | numeric | 0 | 116.13 | 341.25 |
| asp_p50 | numeric | 0 | 174.04 | 351.06 |
| asp_p75 | numeric | 0 | 232.57 | 355.19 |
| asp_p90 | numeric | 0 | 273.54 | 358.60 |
| cCv_p10 | numeric | 0 | 10.84 | 95.00 |
| cCv_p25 | numeric | 0 | 22.00 | 100.00 |
| cCv_p50 | numeric | 0 | 36.69 | 100.00 |
| cCv_p75 | numeric | 0 | 49.42 | 100.00 |
| cCv_p90 | numeric | 0 | 58.46 | 100.00 |
| cHgt_10 | numeric | 0 | 9.95 | 36.00 |
| cHgt_25 | numeric | 0 | 12.13 | 36.00 |
| cHgt_50 | numeric | 0 | 15.35 | 40.00 |
| cHgt_75 | numeric | 0 | 18.31 | 40.00 |
| cHgt_90 | numeric | 0 | 20.04 | 40.00 |
| cat_Cov | integer | 0 | 597 686.44 | 2 999 999.00 |
| dem_p10 | numeric | -11 | 1 071.91 | 3 566.37 |
| dem_p25 | numeric | -7 | 1 118.29 | 3 611.00 |
| dem_p50 | numeric | -6 | 1 183.86 | 3 701.17 |
| dem_p75 | numeric | -6 | 1 253.49 | 3 794.00 |
| dem_p90 | numeric | -6 | 1 308.16 | 3 919.00 |
| dmg_Cov | integer | 0 | 90 426.05 | 1 074 988.00 |
| dmg_wCv | integer | 0 | 66 670.52 | 963 963.00 |
| mnCt_Cv | integer | 0 | 427 530.22 | 2 820 898.00 |
| slp_p10 | numeric | 0 | 10.54 | 41.61 |
| slp_p25 | numeric | 0 | 14.41 | 48.66 |
| slp_p50 | numeric | 0 | 19.01 | 65.91 |
| slp_p75 | numeric | 0 | 23.60 | 73.24 |
| slp_p90 | numeric | 0 | 27.48 | 82.84 |
Italy is divided into regions and each region into provinces. We have a full cover of forest categories of two regions, Veneto and Friuli Venezia Giulia, and the province of Trento. Categories slightly differ in classification, but discrepancies were merged into a single catogory.
Classification can provide us with a binary result: positives/negatives = damaged/undamaged forest areas - with the following alternative outputs per each unit:
Forzieri, Giovanni, Matteo Pecchi, Marco Girardello, Achille Mauri, Marcus Klaus, Christo Nikolov, Marius Rüetschi, et al. 2020. “A spatially explicit database of wind disturbances in European forests over the period 2000 - 2018.” Earth System Science Data 12 (1): 257–76. https://doi.org/10.5194/essd-12-257-2020.